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FOR  DNS/LES  ON  NON-UNIFORM  GRID 
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Dipcirtimento  di  Ingegneria  Aerospaziale 
Seconda  Universita  di  Napoli,  Aversa,  Italia 

Abstract 

This  paper  is  concerned  with  the  development  of  a fourth  order  Finite  Volume 
scheme  for  the  numerical  solution  of  the  incompressible  Navier-Stokes  equations 
on  non-uniform  grids.  In  fact,  the  use  of  non-uniform  computational  grids  is 
inevitable  in  handling  non-homogeneous  flow  computations,  while  numerical 
simulation  of  turbulent  flows  demand  for  higher  order  schemes.  The  effective 
high  order  accuracy  is  obtained  by  reformulating  the  momentum  equation  in 
terms  of  a fourth  order  deconvolved  velocity  field.  Both  a proper  integration  and 
flux  reconstruction  is  implemented  for  the  space  discretization.  A Fractional 
Time-Step  method  for  the  pressure-velocity  de-coupling  is  adopted  and  a second 
order  semi-implicit  scheme  is  used  for  the  time  integration.  Particular  attention 
has  been  devoted  in  developing  congruent  time-accurate  intermediate  boundary 
conditions  for  the  predictor  step. 

1.  Introduction 

Direct  Numerical  Simulations  (DNS)  as  well  as  Large  Eddy  Simulation  (LES) 
demand  for  accurate  and  efficient  numerical  schemes  due  to  the  wide  range  of 
length  scales  involved  in  a turbulent  flow.  In  fact,  low  order  methods  show  high 
numerical  errors  in  the  smallest  resolved  scales.  As  a matter  of  fact,  for  a long 
time,  the  numerical  simulation  of  turbulent  flows  has  been  carried  out  by  means 
of  second  order  Finite  Difference  (FD)  central  scheme  since,  from  the  LES  point 
of  view,  one  performs  an  implicit  application  of  the  top-hat  spatial  filter. 
Moreover,  LES  on  non-uniform  grids  were  initially  performed  in  a 
straightforward  manner  without  taking  into  account  for  the  existing  commutation 
error,  while  only  recently  the  correct  equations  for  LES  on  non-uniform  grids 
were  analysed  [1]  and  much  more  importance  was  given  to  the  correlation 
between  numerical  errors  and  modelling  ones  (e.g.:  [2]).  More  recently, 
conservative  fourth  order  FD  schemes  were  proposed  over  both  staggered  and 
co-located  non-uniform  grids  (e.g.:  [3]). 

As  it  regards  with  the  Finite  Volume  (FV)  method,  the  integral  form  of  the 
Navier-Stokes  equations  appears  the  most  opportune  by  a physical  point  of  view, 
allowing  mass  and  momentum  to  be  a-priori  conserved.  In  this  framework,  an 
evolution  equation  for  the  volume-averaged  field  v is  solved,  obtaining  a second 
order  approximation  for  the  point-wise  velocity  v,  even  by  adopting  higher-order 
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fluxes  integration.  Very  recently,  a fourth  order  FV  compact  scheme  was 
proposed  in  [4],  where  a deconvolution  technique  was  proposed  in  order  to 
recover  fourth-order  accuracy  from  the  volume  averaged  velocity  field.  Actually, 
such  procedure  was  applied  only  in  a post-processing  step  while  the  solved 
variable  remains  the  second  order  averaged  one. 

On  the  contrary,  in  this  paper,  we  follow  the  approach  already  introduced  in  the 
framework  of  the  so-called  Implicit  Structural  Models  [5,  6],  In  that  approach,  a 
modified  integral  equation,  governing  the  evolution  of  an  effective  fourth-order 
variable,  is  obtained  by  means  of  a de-convolution  procedure  applied  on  the 
original  FV  equations.  Following  such  guidelines,  the  development  of  a high 
order  FV  scheme  on  non-uniform  grid  is  illustrated  in  the  framework  of  the 
Fractional  Time-Step  (FTS)  method  for  pressure-velocity  de-coupling.  The 
spatial  discretization  is  performed  according  to  the  Simpson  integration  rule 
while  proper  Lagrange  interpolation  is  used  for  the  fluxes  reconstruction.  The 
time  integration  is  performed  by  means  of  the  semi-implicit  Adams- 
Bashforth/Crank-Nicolson  (AB/CN)  scheme.  Finally,  time-accurate  boundary 
conditions  to  be  associate  to  the  predictor  equation  were  developed  in  a manner 
consistent  with  the  adopted  time  integration.  The  proposed  method  has  been 
validated  in  the  numerical  simulation  of  both  the  two-dimensional  Taylor 
decaying  vortex  solution  and  a time  evolving  mixing  layer. 


2.  Deconvolved  Navier  Stokes  Equations  On  Non-Uniform  Grids 

Consider  the  Navier-Stokes  equations  for  incompressible  isothermal  flows  in  a 
bounded  domain  V,  written  in  integral  non-dimensional  form  over  a Finite 
Volume  (FV)  il(x)cV,  centred  in  x,  whose  boundary  is  denoted  by  3fi(x): 


Jn-  v dS  =0  (1) 

an(x) 


1 


jn¥dS=0 

dft(x) 


(2) 


being  v the  local  volume  averaged  velocity,  l£2(x)l  the  measure  of  the  FV,  n the 
local  unit  vector  outward  to  the  boundary  dQ.  and  F the  momentum  flux  tensor, 
expressed  as  vv  + Ip- Vv/Re.  A proper  initial  field  v0  and  boundary  conditions 


\b  on  dV  must  be  associated  to  the  system  (l)-(2). 


Following  the  formulation  proposed  in  [5],  an  m-th  order  Taylor  expansion  for  v 
around  the  FV  centre  x is  performed,  being  m an  even  integer.  Owing  to  the  cell 
symmetry,  one  gets  (for  sake  of  brevity,  time  dependence  is  omitted): 

v(x)  = (/x-R'm))v  +o(r2)sC;”v  +o{hm+2)  (3) 

being  h a linear  extension  of  the  FV,  e.g.  /z=in(x)ll/3.  The  differential  operator 
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/- th  order  three-dimensional  derivative  in  a Cartesian  reference  system  and 
I 3 

C;w(x)=  II  J {x'k-xkfdx  the  coefficients  of  the  moments  of  the 

PvM  *=i  n(x) 

Taylor  expansion  terms.  By  truncating  and  inverting  Eq.(3),  one  obtains  the  m-th 
order  de-convolved  velocity  v(x)s  [G*m)fv,  whose  properties  in  both  physical 

and  Fourier  space  were  analysed  in  [5].  If  the  inverse  operator  [g^1]  ' is  applied 
on  Eq.(2),  one  gets  an  evolution  equation  for  the  de-filtered  field  v , but 
commutation  terms  appear  for  non-uniform  grids.  Herein,  in  order  to  avoid  the 
explicit  computation  of  such  terms,  the  LHS  of  Eq.(2)  is  simply  re-written,  so 
that: 

= (4) 
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As  a matter  of  fact,  from  the  LES  point  of  view,  Eq.  (4)  should  be  supplied  with 
a suitable  Sub-Grid  Scales  (SGS)  model,  in  order  to  express  the  RHS  in  terms  of 
the  resolved  variable  v . However,  in  this  paper,  we  just  consider  F(v)~F(v) 


without  addressing  this  issue.  In  the  framework  of  Implicit  Structural  Models 
[6],  this  can  be  interpreted  as  an  LES  approach  for  the  top-hat  filtered  variable, 
supplied  by  a sort  of  generalized  scale  similarity  SGS  model.  Finally,  the  de- 
convolution  order  is  fixed  to  m=2  (G2  =G[2))  in  order  to  get  v representing  an 
effective  fourth-order  approximation  to  the  point-wise  velocity  field. 


3.  A Fourth  Order  Deconvolution-based  Scheme  on  Non-Uniform  Grids 

In  this  section,  a fourth  order  FV  method  is  developed  for  2-D  flows  simulation 
on  Cartesian  non-uniform  grids.  The  de-coupling  between  the  velocity  and  the 
pressure  gradient  is  performed  according  to  the  FTS  method  [7],  while  the  semi- 
implicit  Adams-Bashforth/Crank-Nicolson  second  order  scheme  is  adopted  for 
the  time  integration.  Furthermore,  let  us  assume  the  computational  domain 
V = [(kZjxfCfL,]  . Owing  to  its  computational  simplicity,  a co-located 
arrangement  of  the  variables  was  adopted,  hence  the  flux  vectors  defined  onto 
the  face-nodes  must  be  approximated  in  terms  of  the  balanced  variables  in  the 
FV  centre  nodes. 


3.1.  Two-Dimensional  Grid  Definition 

The  grid  points  (see  Fig.  1 ) are  uniformly  distributed  in  x-direction  (supposed  to 
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be  a stream-wise  direction),  i.e.,  x,  ={i-\)hl  + hl/2,  with  h\=L{/Nx  the  step  size, 

i=  1 ,. . .Ni,  being  Ni  the  number  of  FVs  in  x-  direction.  The  y-direction  (supposed 
to  be  a normal  wall  one)  has  non-uniform  grid  spacing  obtained  by  means  of  an 
1 -D  mapping  y = Y {%) , being  If  the  independent  variable  in  the  computational 
domain.  This  latter  is  uniformly  discretized  by  a step  size  H-Lx/Nx,  being  N2  the 
number  of  FVs  in  y-direction.  Furthermore,  for  each  FV  the  face  co-ordinates 
yj=Y(Cj) and  y+j=Y(Cj+l)  , are  defined  for  j=  1...  N2  and  thus,  the  FV  grid  node 

results  y j = (yj  + yj )/  2 , being  yT  = y~_,  for  grid  construction.  The  mesh  size 
in  y-direction  is  defined  as  h^j)^  yj  - y~  = Y(fj+1)-Y(fj)  having  assumed  a 
smooth  mapping  ( h2/H=0(\)  ) so  that,  the  FV  definition  is 
Q..  =[x, -h^/2;xi  + h]/2]x[yj-h2(j)/2-,yj  + h2(j)/2\.  It  is  noteworthy  that  one 
can  simply  express  the  face  co-ordinates  as  xf  = x,  ±/i,  /2,  yj  = y . ± h2(j)/2  . 

3.2.  The  FTS  Procedure  and  The  Discrete  Time  Integration 
In  Eqs.(l)  and  (4)  the  diffusive  terms  along  the  y-direction  are  integrated  in  time 
by  means  of  the  Crank-Nicolson  scheme,  while  the  Adams-Bashforth  one  is 
adopted  for  all  the  other  terms.  In  the  present  FTS  method  {pressure-free 
projection  method ),  first  an  equation  for  a non-solenoidal  vector  v*  is  obtained 
by  integrating  Eq.(4)  and  eliminating  the  pressure  term: 


In  the  previous  relations,  the  operator  D was  split  as  D = Dx+  Dy,  along  the 
Cartesian  directions,  being  x1  and  y±  the  face  co-ordinates  of  Q.  Observe  that  the 
de-convolution  procedure  does  not  increase  the  computational  cost  since  a semi- 
implicit  procedure  has  been  adopted  for  the  time  integration.  Therefore,  the  LHS 
will  simultaneously  take  into  account  for  both  deconvolution  and  time 
integration 

The  predicted  velocity  field  v*  must  be  corrected  by  means  of  a pure  gradient 
field  according  to: 

vn+1  =\*  -V0  . (7) 

Therefore,  once  the  vector  field  v*  was  computed,  the  continuity  constraint  is 
enforced  at  the  new  time  level  f+I  by  means  of  the  projection  step: 
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; ^ = n.(v--vr')onav  ,8) 

associated  to  the  proper  normal  boundary  conditions. 

3.3.  Fourth  Order  Spatial  Discretization 

The  achievement  of  a real  fourth  order  space  accuracy  is  obtained  by  means  of 
the  Simpson  integral  discretization  along  with  explicit  Lagrangian  interpolation 
for  the  fluxes  discretization.  Although,  at  the  same  accuracy,  a Lagrangian 
polynomial  involves  a wider  stencil  than  a Hermitian  interpolation  (as  recently 
proposed  in  [4]),  the  former  approach  remains  simply  applicable  in  the  non- 


uniform  direction. 


operator 


g2  = /x+- 


h2Ay ) a2 
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discretized  to  fourth  order  accuracy  by  considering  second  order  central 
difference  formulas  for  the  spatial  derivatives.  Having  deconvolved  the  velocity 
field  to  fourth  order  accuracy,  the  integral  fluxes  in  Eqs.(6)  and  (8)  are 
congruently  discretized  by  means  of  the  Simpson  formula.  In  such  a way,  one 
has,  as  an  example,  for  the  net  flux  along  x-direction: 

y-j  9 ' i.i 

wherein,  the  face-node  unknowns  (see  Fig.l)  are  expressed  in  terms  of  the  grid 
nodes  values.  A high-order  Lagrangian  interpolation  procedure  is  adopted  by 
factorising  the  function  along  each  direction  and  approximating  both  the  factors 
by  means  of  two  third  degree  polynomials,  i.e.:  f(x,y)  = f{x,y)  = Ll(x)L2(y). 
This  way,  one  obtains  a global  25  grid-nodes  computational  molecule. 


3.4.  Boundary  Conditions 

The  original  system  (1),  (4)  was  de-coupled  in  the  separate  prediction  (6)  and 
projection  (8)  equations.  The  solution  of  this  latter,  together  with  Eq.  (7),  allows 
the  intermediate  velocity  v*  to  be  corrected  by  imposing  the  exact  normal 
velocity  component  n-v;/'+1  on  dV.  Since  the  projection  can  not  correct  the 
tangential  velocity  component,  this  latter  must  be  congruently  assigned  for 
solving  Eq.(6).  In  this  paper  we  propose  a procedure  for  assigning  time  accurate 
boundary  conditions.  In  fact,  by  taking  the  limit  of  Eq.(6)  for  vanishing  grid 
spacing  ( G2—>I  for  h—>0)  and  projecting  it  along  the  direction  tangential  to  the 
boundary,  one  gets: 

JV"+AL^L)  + (10) 

r 2Re  3y2  j {'  2Re  L 


-tjY-[(3vvr-(wr']-^[3 
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When  properly  spatially  discretized,  Eq.(10)  provides  the  correct  second  order 
time  accurate  solution  on  the  frontier,  as  illustrated  in  the  next  section. 

4.  Numerical  Results 

The  adopted  test  case  is  the  classical  2-D  Taylor  decaying  vortex  solution  at 
Re=l.  The  error  estimations  are  performed  in  the  computational  domain  V=[-7t, 
ri\Y.[-7Z,  tc_ |,  by  taking  the  L„  norm  of  the  difference  between  the  x-component  of 
the  exact  and  numerical  velocity  field.  In  a first  test,  for  bi-periodical  boundary 
conditions,  the  effect  of  the  de-averaging  procedure  was  studied,  having  used  the 
fourth-order  flux  integration  (9).  The  resulting  errors,  computed  after  100  time 
steps  (At^lO-4),  are  shown  in  Fig.2  versus  the  normalised  grid  size  in  a double 
logarithmic  scale,  for  both  uniform  and  non-uniform  grids.  It  can  be  noted  that 
only  in  the  presence  of  the  de-averaging  procedure  (m= 2)  an  effective  fourth 
order  accuracy  is  reached  in  computing  unsteady  solutions.  The  space  accuracy 
is  also  checked  for  Dirichelet  boundary  conditions  in  the  y-direction,  as 
illustrated  in  Fig.3  on  both  uniform  and  non-uniform  grids.  Moreover,  the 
correctness  of  the  proposed  boundary  conditions  (10)  is  analysed  from  the 
results  in  Fig.4,  where  the  time  accuracy  tests  is  reported.  A single  time 
integration  was  conducted  in  order  to  avoid  stability  problems  when  working 
with  high  time  steps.  The  errors  are  reported  versus  the  used  time  step,  showing 
a third  order  slope  according  to  the  fact  we  are  evaluating  the  direct  error  on  a 
single  time  step,  that  corresponds  to  a second-order  local  truncation  error. 

Finally,  the  procedure  was  tested  for  a time  evolving  mixing  layer.  The  initial 
configuration  is  the  same  adopted  in  [9],  i.e.  a hyperbolic  tangential  velocity 
profile  u(y)  = tanh  2y/Sl  (being  2 Sx  the  initial  vorticity  thickness) 
submitted  to  a white  noise  perturbation  plus  a deterministic  sine  perturbation  at 
k4=(2n/X:i)  with  A,a=7  Sx . Kelvin-Helmholtz  instabilities  lead  to  the  development 
of  vortices  which  in  a later  stage  roll-up  and  merge.  This  is  a good  example  for 
the  tendency  of  2D  turbulence  to  transfer  energy  from  small  to  large  scales  thus 
requiring  an  accurate  numerical  simulation.  The  preliminary  results  of  this  study 
are  reported  in  Fig.s  5 for  a 1282  computational  grid  on  a domain  V=[0,  4A.a]x[- 
2X,a,  2A,a,]  at  Re  = /v  = 250  showing  the  salient  features  of  backscatter 
transfer  energy  and  the  appearance  of  a kA  range  according  to  the  LES  performed 
in  [9]. 
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Figure  I : description  of  the  adopted  2-D  grid, 
non-uniform  in  the  y-  direction.  For  each  j-th 
finite  volume,  the  vertical  flux  section  is 
centred  with  respect  to  the  node  j i.e., 
.y±(yj)=yJ±h2(j)l2-  The  unknown  variables  u,  v, 
p are  co-located  on  the  grid  nodes  (•). 


Figure  2:  Space  accuracy  tests  with  Bi- 
periodic  boundary  conditions  applied.  Effects 
of  the  de-averaging  procedure  on  the  accuracy. 
Errors  computed  after  100  time  steps  on  both 
uniform  and  non-uniform  grids. 


Figure  3:  Space  accuracy  tests  with  Periodic-  Figure  4\  Time  accuracy  tests.  Errors 
Dirichlet  boundary  conditions  applied.  Errors  computed  after  one  time  steps  on  both  uniform 
computed  after  100  time  steps  on  both  and  non-uniform  grids  of  602  CV’s. 
uniform  and  non-uniform  grids. 
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